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Phase separation in salt-free symmetric mixtures of oppositely charged rodlike 
polyelectrolytes is studied using quasi-analytical calculations. Stability analyses for 

the isotropic-isotropic and the isotropic-nematic phase transitions in the mixtures are 

CM- 

,_h ' carried out and demonstrate that electrostatic interactions favor nematic ordering. 

Coexistence curves for the symmetric mixtures are also constructed and are used 

00 : 

t— 1 , to examine the effects of linear charge density and electrostatic interaction strength 

on rodlike polyelectrolyte complexation. It is found that the counterions are uni- 
formly distributed in the coexisting phases for low electrostatic interaction strengths 
dictated by the linear charge density of the polyelectrolytes and Bjerrum's length. 
However, the counterions also partition along with the rodlike polyelectrolytes with 
an increase in the electrostatic interaction strength. It is shown that the number 

O : 

density of the counterions is higher in the concentrated (or "coacervate" ) phase than 
T — 1 . in the dilute (or supernatant) phase. In contrast to such rodlike mixtures, flexible 

>: 

polyelectrolyte mixtures can undergo only isotropic-isotropic phase separation. A 

^" : 

comparison of the coexistence curves for weakly-charged rodlike mixtures with those 

^ '. 

^. . of analogous flexible polyelectrolyte mixtures reveals that the electrostatic driving 

o 

force for the isotropic- isotropic phase separation is stronger in the flexible mixtures. 

^H . 
>'■ 

[Xj ■ I. INTRODUCTION 

For numerous biological processes^- and emerging technologies^, complexation between 
oppositely charged polyelectrolytes 7-13 is the underlying fundamental phenomenon. How- 
ever, our understanding of how the electrostatic attraction between opposite charges on the 
polyelectrolytes is coupled with the conformational characteristics of polyelectrolyte chains 
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is very limited. Motivated by a plethora of relevant biological processes and the impetus 
for developing advanced technologies such as underwater adhesives 5 and biosensors^, exten- 
sive experimental^ and theoretical^— work has been carried out to understand the 
mechanisms of polyelectrolyte charge complexation. 

A few insights into the origin of polyelectrolyte complexation have been provided by sim- 
ulations. Langevin dynamics simulations 1 - of two oppositely charged polyelectrolyte chains 
have revealed that the entropy of counterion release is the predominant driving force for 
complexation in highly charged polyelectrolytes. In contrast, the simulations show that 
direct electrostatic attractions drive the association of two weakly charged polyelectrolytes 
whose counterions are not condensed onto the chains. Such simulations involving only a 
single pair of polymers are of course relevant only to extremely dilute solutions of polyelec- 
trolytes and are only suggestive of complexation phenomena at finite concentrations. At 
such higher concentrations, complexation between oppositely charged polyelectrolytes has 
the character of a phase separation process in which a supernatant phase that is extremely 
dilute in polymer macroscopically separates from a polyelectrolyte rich phase that is either a 
fluid ("complex coacervate") or a solid ("precipitate"). In this paper we shall be exclusively 
concerned with systems that produce fluid-like polyelectrolyte complexes properly classified 
as complex coacervates.-^&i^r— Extending conventional Langevin, Monte-Carlo, or molec- 
ular dynamics particle-based simulations to polyelectrolyte mixtures at finite concentrations 
has been hindered by the need for extremely large computational resources to address the 
long-range character of the electrostatic interactions and the inherently slow kinetics of the 
dense coacervate phase and phase separation process. 

Development of a theory to study the phase separation in a mixture of oppositely charged 
polyelectrolytes is a daunting problem due to the non-trivial coupling between long range 
effects such as electrostatics and chain connectivity, and short range excluded volume inter- 
actions. The theoretical treatment is further complicated by the inevitable need to include 
electrostatic screening and correlation effects in the theory, analogous to the classical Debye- 
Hiickel theory of electrolyte solutions, to capture the electrostatic forces that drive phase 
separation. This is the reason behind the inability of the self-consistent field theory (SCFT), 
a type of mean-field theory, to capture even the signature of a phase separation^— . In other 
words, in a field-theoretic description of a mixed polyelectrolyte system, one needs to go be- 
yond the mean-field or "saddle-point" approximation in order to study phase separation 



phenomena. Typically, the random phase approximation^!"— is used to compute the leading 
correction to the saddle-point results. In the case of flexible polyelectrolytes^, it can be 
shown rigorously that the random phase approximation is valid in the dual limit of high 
monomer densities and low small ion densities. Physically, such a situation is realized in 
concentrated solutions of weakly charged polyelectrolytes. More sophisticated field-theoretic 
results beyond the random phase approximation have also been obtained (including a phase 
diagram) for a particular model of flexible polyelectrolyte mixtures using the self-consistent 
one-loop (Hartree) approximation 24 . Most recently, direct numerical simulations of a re- 
lated field theory model (so-called "field-theoretic simulations") have been carried out^ 1 ^ 
to capture the effects of fluctuations and correlations to all orders, without invoking any 
approximation. An important result that came out of these simulations of flexible poly- 
electrolyte mixtures is that the phase boundaries in the concentrated regime are accurately 
predicted by the random phase approximation. 

In this work, we consider solutions containing binary mixtures of oppositely charged rod- 
like polyelectrolytes and study their phase behavior using the random phase approximation, 
which, as noted above, has been previously shown to provide an accurate description of 
complexation phenomena in flexible polyelectrolyte mixtures. An important focus of our 
study is a comparison of the phase behavior of rodlike and flexible polyelectrolyte mixtures. 
The rodlike system is particularly relevant to a recently developed biosensor technology 6 
involving complexation of cationic conjugated polyelectrolytes with anionic DNA. 

A fundamental question that we have tried to answer is: what is the role played by the 
flexibility of the polyelectrolytes in the phase separation processes that lead to a complex 
coacervate? Unlike flexible polyelectrolyte mixtures^"—, local enhancements in concentra- 
tion of rodlike systems can produce orientationally ordered liquid crystalline phases (such 
as nematic, smectic and cholesteric phases).— Thus, a coacervate produced by complexing 
oppositely charged rodlike polymers may be a liquid crystal rather than an isotropic fluid. 
In this paper, where we consider "symmetric" achiral rodlike polymers with fore-aft symme- 
try, equal length, and equal but opposite charge, we have explored the competition between 
isotropic and nematic ordering in rodlike polyelectrolyte mixtures and elucidated how the 
phase boundaries depend on different parameters of the system. 

Before presenting our theory and numerical results on the binary mixtures, it is worthwhile 
to briefly review the theoretical literature aimed at understanding the phase separation of 



single component flexible^— and rodlike^ - — polyelectrolyte solutions. These studies have 
been primarily focused on poor solvent conditions. In the case of solutions containing a 
flexible polyelectrolyte, it has been shown that the dependence of the phase boundaries 
on the polyelectrolyte chain length is weak (referred to as the "weak polymer effect"—). 
Note that the condition of Donnan equilibrium was invoked in these calculations so that 
locally the system is electroneutral. This implies that a local increase in the charge on the 
polyelectrolyte leads to an increase in the local concentration of counterions. The localisation 
of the counterions is entropically unfavorable and hence, the region of two-phase coexistence 
shrinks with an increase in the charge on the polyelectrolyte chains^"— . 

Phase-separating mixtures of two oppositely charged poly electrolytes are fundamentally 
different from these solutions containing a single kind of polyelectrolyte in two ways. Firstly, 
oppositely charged polyelectrolyte mixtures can phase separate even under good solvent 
conditions for the two individual polymer components. Secondly, it is clear in the binary 
mixture case that counterions are free to partition very differently upon phase separation. In 
particular, Donnan equilibrium can be satisfied by bringing the oppositely charged macro- 
molecules together in the concentrated coacervate phase,— - — while freeing the small ions 
to gain entropy by populating both dilute and concentrated phase in charge compensating 
proportions. Clearly, it is necessary to account for these extra degrees of freedom available 
to the counterions in studying phase separation in mixtures of oppositely charged polyelec- 
trolytes. 

In the case of single component rodlike polyelectrolyte solutions, it has been shown 
recently^— that electrostatic interactions among like-charged rods favors orientational or- 
dering at concentrations lower than the overlap concentration. This non-trivial result is 
an outcome of long-range, multi-rod correlations present in these systems. A consequence 
is that, in addition to the lyotropic isotropic-nematic transition observed in neutral rod- 
like solutions, thermotropic isotropic-nematic as well as nematic-nematic phase transitions 
have been predicted for charged rodlike polyelectrolytes. Note that these predictions are 
a consequence of collective phenomena at higher concentrations that are not manifest in 
earlier theories of dilute rodlike polyelectrolytes, where the twisting effect^ 1 ^ of pairwise 
electrostatic interactions between rods was shown to destroy nematic ordering^"—. 

With this background on the phase separation in polyelectrolyte solutions, we present 
our theoretical model for the binary mixtures. This paper is organized as follows: the theory 



is presented in Section HI1 the results are presented in Section IHIl and Section [IV] contains 
our conclusions. 

II. THEORY 

We consider a symmetric binary mixture of rodlike polyelectrolytes bearing opposite 
charges on their backbone in the presence of the counterions originating from the poly- 
electrolyte chains. By symmetric mixtures, we mean mixtures containing equal number 
of polyelectrolyte chains, which are identical to each other in every aspect other than the 
sign of the charge they are bearing. For generality in the following, we maintain notation 
representing different polymeric species such as n + and n_ representing the number of poly- 
cationic and polyanionic chains, resepectively. For the theoretical treatment, we consider a 
binary mixture of rods with the same diameters (= d) for each species but lengths L& for 
k = +, — . Furthermore, er^ is taken to be the linear charge density for polyelectrolytes of 
type k = +, — so that there are n c k = a^Lk counterions released by each polyelectrolyte of 
type k. As the charge on the polymers is spread uniformly along their length, this descrip- 
tion is often referred to as the "smeared charge" model. We specialize to the symmetric 
case only as needed after plugging n + = n_ = n, L + = L_ = L and a + = o~_ = a at the 
appropriate place. Also, we work in the canonical ensemble with the volume of the system 
denoted by Q. 

The solvent is treated implicitly and as a uniform dielectric background for the purpose 
of computing electrostatic interactions. The small ions are taken to be point like and we 
ignore short range structure, correlations, and polarization effects in the solvent. 

A. Qualitative Picture and Scaling Analysis 

Before presenting quantitative details of our theory for the symmetric rodlike polyelec- 
trolyte mixtures, a qualitative picture can be drawn by using three important results avail- 
able from the literature. First, the macrophase separation in the isotropic phase (by anal- 
ogy with symmetric mixtures of flexible polylectrolytes) takes place at very low monomer 
densities^—, when the electrostatic attractions between oppositely charged polyelectrolytes 
are strong enough to drive the phase separation. Second, the isotropic-nematic phase transi- 
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tion observed in neutral rodlike molecules and treated by Onsager and many others^^^r— 
takes place at relatively high monomer densities above the overlap concentration and is 
driven primarily by entropic effects associated with the excluded volume of rodlike objects. 
Thirdly, we recall that addition of charges to rodlike polymers can induce thermotropic 
(temperature dependent) isotropic-nematic phase transitions, as opposed to the lyotropic 
(concentration dependent) transitions familiar in neutral rodlike polymer systems. The ori- 
gin of this thermotropic behavior lies in the long-ranged electrostatic interactions, which 
depend on the charge on the polyelectrolytes and the temperature dependent Bjerrum's 
length Ib = e 2 /(47T€o€ r kBT), which is the distance at which the Coulomb interaction energy 
between two elementary charges is of the order of thermal energy (= ksT, ks being the 
Boltzmann's constant and T is the temperature). Here, eo and e r are the permittivities of 
vacuum and the medium, respectively. 

By integrating these three competing effects, a qualitative phase diagram of a symmet- 
ric rodlike polyelectrolyte mixture is sketched in Fig. Q] in the coordinates of Ib/1 versus 
2nL 2 d/Q, where I is a reference monomer length, L is the overall rod length, d is the rod 
diameter, 2n is the total number of cationic and anionic rods, and Q is the system volume. 
Broadly, we see two trends in the figure. At high temperatures, i.e. low values of Ib/1, 
non-interacting isotropic-isotropic and isotropic-nematic coexistence regions appear at low 
and high rod concentrations, respectively. As the temperature is lowered {Ib/1 raised), these 
features collide and ultimately at low temperature there is a broad region of two phase 
coexistence between a dilute isotropic supernatant phase and a semidilute or concentrated 
nematic "coacervate" phase. At intermediate values of Ib/1, more complex phase behavior 
is evident, including the possibility of three-phase coexistence between two isotropic phases 
and a nematic phase. 

With regard to the isotropic-isotropic coexistence region in Fig. [TJ which is also present 
in flexible polyelectrolyte mixtures, a scaling analysis is helpful to understand the variables 
that control its extent and shape and differentiate the cases of rigid and flexible systems. For 
this purpose, we begin by considering a symmetric mixture containing an equal number n 
of oppositely charged flexible polyelectrolytes in the absence of any counterions. Physically, 
such a situation is realized in a solution containing polyacids and polybases. Apart from 
electrostatic interactions, excluded volume effects can be simply accounted for (in an implicit 
solvent model) by an excluded volume parameter w. This particular system was recently 



studied using the field-theoretic simulation technique 2 ^ 1 ^. The osmotic pressure of such a 
system in the dilute regime 2 ^^, i.e. p = 2nN/Q <C p*j> ~ i/(47tZb<7 2 ./V 2 ), with p, p*? being 
the monomer density and the overlap concentration, respectively, is given by 
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where K 2 d = AkIb^N p is the Debye screening length in the asymptotic dilute regime where 
the polyelectrolytes behave as multivalent macroions. We note that although the polymers 
are assumed to be flexible, the expression that we use for pj assumes that they adopt 
extended rodlike conformations at infinite dilution. In the above expression, a is the charge 
per monomer (in units of the fundamental charge e), which is the same for each kind of 
polyelectrolyte in the assumed symmetric mixture and N is the total number of monomers. 
In Eq. (U the first and second terms correspond to the translational entropy of the 
polyelectrolyte chains and the effect of excluded volume interactions on the pressure, re- 
spectively. The third term containing Kd, a Debye-Hiickel contribution, originates from 
attractive electrostatic correlations among dilute, non-overlapping polyelectrolyte chains. 
In contrast, at higher polymer concentrations where chains overlap the internal structure of 
the polyelectrolytes plays a significant role. The dominant electrostatic contribution to the 
osmotic pressure in the high density, overlapping regime (p ^> p**) can be represented by the 
scaling expression 

pe 3 
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where F{x) is a dimensionless scaling function that we assume to be a power law ~ x m for 
large argument. The exponent m is determined by the requirement that in the concentrated 
regime, the local properties are not different for a solution having multiple chains containing 
N monomers each or a single chain that fills space with the same overall segment density p 
(N — > oo). In other words, the electrostatic contribution to the osmotic pressure must be 
independent of iV at fixed p, i.e. Pf ~ iV°. This requirement, along with the fact that the 
overlap concentration scales as pj ~ l/(47r/scr 2 A r2 ), leads to the well-known result^ 
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where k 2 = Airls^p. In the concentrated regime, the osmotic pressure of flexible, symmetric 
polyelectrolyte mixtures is thus dictated by a balance between excluded volume and this 



modified electrostatic correlation energy, i.e. 
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We can now retrace these scaling arguments for the case of symmetric mixtures of rod- 
like polylectrolytes. The osmotic pressure in the dilute concentration regime (p <C p* ~ 
2nL 2 d/Q) is unchanged except that we replace the excluded volume coefficient w with the 
Onsager expression irdl 2 /2, where / = L/N is a monomer segment length, 

/I r ri/1„ V / 
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In the concentrated regime, the electrostatic free energy of a system of rods can be ap- 
proximated by the sum of the self-energies of individual rods experiencing the electrostatic 
potential induced by the other rods. For the purpose of estimation, these are placed on a 
periodic lattice. It is well-known that the self-energy of a charged cylindrical rod of length 
L and linear charge density a is divergent^^, and is given by 

U cyl r 2 - -^ 
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where 5 is a cut-off introduced to regularize the self-energy and R is the radius of a unit cell 
(Wigner-Seitz) applied to each cylinder. Using the fact that R ~ 1/y 7 /? in the concentrated 
regime and writing the electrostatic free energy of the solution containing n rods as F e = 
nlI C yi, the electrostatic contribution to the pressure is given by 
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Comparing Eqs. |3] and [TJ it is clear that the qualitative difference between flexible 
and rodlike polyelectrolyte mixtures is in the functional form for the attractive electrostatic 
contribution to the osmotic pressure. In the case of rodlike polyelectrolytes, this term scales 
like Kp ~ p in contrast to /%/ ~ p 3 / 4 for flexible coils. Thus, at the same set of electrostatic 
parameters and monomer densities, and for weakly charged polyelectrolytes where k p < 1, 
electrostatic correlation effects are weaker for rodlike polyelectrolytes in comparison with 
flexible polyelectrolyte mixtures. Since the electrostatic contribution to the pressure is 
responsible for polyelectrolyte complexation, it follows that rodlike polyelectrolyte mixtures 
are less prone to complex coacervation than analogous flexible polyelectrolyte mixtures. 
The regions of two-phase coexistence sketched in Fig. [TJ should thus be narrower in rigid 
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rod systems. This general observation will be born out by detailed numerical calculations 
in Section IIHI of the spinodals and binodals of each type of system based on the free energy 
expression presented below. 

B. Free Energy 



An expression for the Helmholtz free energy of the assumed mixture of rodlike polyelec- 
trolytes, counterions, and implicit solvent is derived in Appendix A using Onsager's treat- 
ment for the neutral interactions and the random-phase approximation for the electrostatic 
interactions. Explicitly, the free energy is given by 
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Here, k is the inverse of the Debye screening length for small ions given by k 2 = 
4tt^s 2~2-y=c+ c- Z^n^/Vt so that Z 7 is the valency (with sign) of the charged species of type 
7. We define a similar object for the two polymer species: n 2 . = 4iil B ZlalLln k /Q, Z k be- 
ing the valency of the charged monomers of type k. The functions p k (vL,n) represent the 
probability distribution function for finding a rod of type k oriented along the unit vector u 
when the director is chosen to be the unit vector n. Each distribution function satisfies the 
normalization condition j dup k (u, n) = 1. Finally, the primed superscript in the expression 
for F w indicates that the j = k terms are omitted from the double sum. 

So far, the theory is quite general and the subsequent analysis depends on the functional 
form for the orientational probability distribution function for the rods. In this work, we 
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focus on possibilities for phase separation that involve the isotropic and nematic phases only. 
For these two phases, p k (u,ii) = p k (u ■ n) and the free energy for the each phase can be 
evaluated as described below. 

C. Isotropic Phase 



We begin by considering the isotropic phase, in which case the orientation distribution 
function for the rods is independent of the angle between the director n and the unit vector 
along the axes of the rods (u). In this case, pk(u-h) = I/Att (obtained from the normalization 
condition) for k = +, — so that the free energy in Eq. [S]can be written as 
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- M is the sine integral. 



where Si(x) = L dt- 1 

In the following, we rigorously show that the scaling argument presented in section III Al 
predicting a logarithmic correction to the free energy of isotropic phase in the mixtures of 
rods is indeed correct. For the analysis, we use the following asymptotic expressions for 
t is °(qL k ) : 

1, qL k ->■ 

-n/qL kl qL k -)• oo. 

Also, we consider two limiting cases of short and long rods. For these cases, the electrostatic 
contribution to the free energy depends on how these limits are being approached. These 
limits can be approached by either fixing the linear charge density or the net charge per 
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rod during the variation of length L k . In the former case, Kk increases (decreases) with an 
increase (decrease) in Lk and diverges (vanishes) strictly in the limit of L k — > oo (L& —¥ 0). 
On the other hand, for the latter approach, Kk is held fixed, while approaching the limits 
of Lfc — > oo or Lk — > 0. Consider the latter case so that Kk is well-defined in approaching 
either limits of L k - 



Defining a Debye-like parameter proportional to the total ionic strength, K 2 ff 



K 2 + 



Ylk=+ - K 1 anc ^ approaching the limit of short rods, L k — > while keeping the charge per 
rod fixed, the free energy becomes 
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Similarly, consider the limit of long rods, i.e. Lk being large and approaching oo. Approach- 
ing this limit while keeping the charge per rod fixed, the free energy becomes 
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where the electrostatic contribution to the free energy is given by 
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Here, we have introduced k 2 = I j^- + -j=- ) . Furthermore, note that in writing Eq. [5TJ we 
have used the aysmptotic form for t tso in the limit of qL — > oo for the entire range of q (e.g., 
even in the case of q = 0). This causes the integral in Eq. [21] to diverge, while the original 
integral in Eq. [16] is convergent. These divergences are mere artifacts of the approximation 
scheme. Despite these divergences, we show that the leading contribution to the free energy 
of the rods is of the form — na 2 l B L\nn/Q, as described using the scaling arguments. 

For weakly charged polyelectrolytes so that ^ < 1, Eq. [2T] can be written as (see 
Appendix B) 
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where the leading term in the free energy expression has the functional form similar to the 
one described in Eq. EJ On the other hand, for strongly charged rods in weakly screened 
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solutions h 2 ^> 1 and k — > 0. In this limit, F e iSO becomes (see Appendix B) 
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Like the weakly charged rods, the free energy of the strongly charged system has the same 
functional form as in Eq. [6j 

Using these approximate expressions for the free energy, the osmotic pressure can be 
readily computed using the thermodynamic relation P = — (dF lso /dQ) nT . For the limiting 
case of L k — > 0, this gives 

k B T 2^ Nk + 2^ P^ 24tt' [ ] 
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which is the well-known Debye-Huckel limiting law for electrolyte solutions. Here, p k = 
n k N k /Q is the monomer number density of type k = +, — and p 1 = n y /Q is the number 
density of counterions of type 7 = c+,c—. 

For the other limiting cases presented above, i.e., for Kp< 1 

P {L k -»■ OO} ^ Pk ^ K 3 7td 2 «J 7T«J 

where 4 = L k /N k is the monomer length. Similarly, for k 2 ^> 1 and k — ¥ 
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Note that the leading contribution to the pressure coming from the electrostatic correlations 
of the charged rods is negative and of the form — k 2 as already anticipated from the simple 
scaling arguments in section III At This particular contribution can drive macrophase sepa- 
ration in isotropic solutions of rodlike polyelectrolytes, a theme that is considered later in the 
paper. Next, we turn to consider the weakly ordered nematic phase using the approximate 
free energy given in Eq. [HJ 
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D. Weakly Ordered Nematic Phase 



For a nematic phase, the distribution function can be reduced to a function of the cosine of 
the angle between the director and rod orientation, i.e. Pk{u, n) = pk{u ■ n) for k = +, — . In 
order to study a fully formed nematic phase, the complete functional form for the probability 
distribution function is needed. However, a stability analysis for a weakly ordered nematic 
phase can be carried out without knowing the probability distribution function a priori by 
assuming a two-term Legendre expansion for Pfc(u • n) written as 
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Here, Sk is the nematic order parameter, given by 
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Here, we have used the normalization condition J dupk(u ■ n) = 1 to write the Legendre 
expansion. Note that Sk = for the isotropic phase and 1 for the completely ordered nematic 
phase, for which pk(u ■ n) = S(u — n). A similar analysis for the case of single component 
rodlike polyelectrolyte solutions has been carried out in Refs.— _ — . 

Using Eq. [571 tk(qLk) for a weakly ordered nematic phase can be written as t k (qL k ) = 
f s °(qL k ) + St k (qL k ), where 5t k (qL k ) is given by 
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Using these expressions for pk and tk{qLk), the free energy of a weakly ordered nematic 
phase (expanded to second order in Sk) can be written as 
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Here, we have used the approximation sin# ~ [1 — 5P2(cos#)/8] vr/4, P2 being the Legendre 
polynomial of order 2, in writing AF W , and 
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Note that the contribution to the free energy coming from electrostatic correlations (i.e., 
AF e ) is negative. In other words, electrostatic correlations favor the nematic phase. This 
observation will prove important to the the stability of the weakly ordered nematic phase 
with respect to the isotropic phase. 

Some useful results can be inferred from the above expression for the free energy of a 
weakly ordered nematic phase. The spinodal of the isotropic-nematic transition (i.e. stability 
limit of the isotropic phase in the nematic region) is given by the condition 
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For the case of symmetric mixtures, i.e. equal length 
density of the rods (= n/fl), Eq. I3T1 reduces to 
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and c p = 2nL 2 d/Q is the dimensionless parameter of the order of overlap concentration of the 
rods. For numerical purposes, dimensionless parameters are introduced after renormalizing 
different parameters with the length of the rods so that k 2 = k 2 L 2 = 47rlB<Jc p Ll/d, k 2 m = 
L 2 {k 2 + + k 2 _) = 47rlBO' 2 CpL 2 l/d, and q = qL. Furthermore, in order to make a qualitative 
connection with the Manning's theory of counterion condensation^ later in this work, we 
have renormalized Bjerrum's length using monomer length (/ = L/N) and defined Ib = Ib/1- 
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Note that for the symmetric mixtures, the spinodal can also be obtained after putting 
S + = £_ = S and then evaluating d 2 F/dS 2 = 0. Furthermore, Sin < corresponds to the 
regime where the isotropic phase is unstable relative to the weakly ordered nematic phase. 
It is clear from Eq. [38] that the electrostatics and the steric effects (last two terms in square 
brackets) act against the orientational entropy to drive the system from the isotropic phase 
to the anisotropic nematic phase. In the absence of electrostatics^ 1 ^, the nematic phase 
becomes stable when c p > 16/tt ~ 5.1. It is also clear from Eq. [3B]that the nematic phase 
in polyelectrolyte mixtures becomes stable at lower polymer concentrations compared to the 
neutral mixtures. 

Numerical solutions of Eq. [38] are presented in Fig. [2] for symmetric mixtures without 
and with the counterions (Fig. |2^l and [2b, respectively) and for different linear charge 
densities. On the right hand side of these boundaries, the nematic phase is stable. It is clear 
from Fig. [2] that the nematic phase becomes stable at lower monomer densities with an 
increase in linear charge density at fixed Bjerrum length. Furthermore, comparing Fig. [2^ 
and [2b, it is evident that the nematic phase becomes stable at higher monomer densities 
in the presence of counterions, i.e. the stable regime of the nematic phase is smaller in the 
presence of counterions. Note that these numerical results are in agreement with earlier 
theoretical predictions for one-component rodlike polyelectrolyte solutions that electrostatic 
interactions favor uniaxial ordering of the rods.— ~— 

E. Nematic Phase for Symmetric Mixtures: Variational Treatment 

In order to study a nematic phase with an arbitrary magnitude of the orientational order 
parameter, we need to resolve the full probability distribution function. In principle, this can 
be done by a calculus of variations approach by minimizing the free energy of the system, 
which leads to an integral equation^. For the case of neutral polymers, the problem has 
been attacked by three different routes. The first is a variational treatment using Onsager's 
trial function^i^^, where the variational parameter is determined by minimizing the free 
energy. A second route is through the use of a Legendre expansionist and determining 
the coefficients in the Legendre series that minimize the free energy. A third scheme is 
to directly attack the integral equation in real spaced using a non-linear equation solver. 
From a computational point of view, the last two routes are more demanding and become 



16 

especially difficult for strongly ordered nematic phases (S > 0.9). The first route is the 
easiest and readily describes a nematic phase with arbitrary order. However, it leads to a 
slight overprediction^i^ of the coexisting densities at the isotropic-nematic transition. 

In this work, we adopt Onsager's variational approach to study nematic phases in the 
symmetric polyelectrolyte mixtures. Due to the assumed symmetry, the probability distri- 
bution functions for the two types of rods must be the same and only one of the distribution 
functions needs to be considered. In Onsager's approach, the probability distribution func- 
tion is taken to be of the form 

a cosh [a (u • n)] 
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(40) 
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where a is a variational parameter, which is determined by minimization of the free energy. 
The order parameter corresponding to this distribution function for a symmetric mixture is 
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so that S = for a = corresponds to the isotropic phase and 5 = 1 when a — > oo, 
corresponding to a perfectly ordered nematic phase. The free energy density can be written 
in terms of a so that the free energy density of the symmetric mixtures becomes / = F/Q = 
fen + f w + fe, where 
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where c s = aLc p is the dimensionless number density of small ions (both positive and 
negative) and I2 is the modified Bessel function of order 2. The integrals in the expression for 
t(q, x) can be readily evaluated using Gauss-Legendre and Gauss- Chebyshev quadratures for 
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x' and y integrals, respectively. Also, the integral over q ranging from to oo can be evaluated 
using Gauss- Legendre quadrature after using the transformation q — (1 — z)/(l + z). 
The free energy density can be optimized with respect to a so that -J- = 0, which gives 
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Note that for a root (= a*) of Eq. H^l to be a minimum of the free energy, Q^\ a = a * > 
must be satisfied. In other words, some of the roots of Eq. SS] may correspond to a local 
maximum in the free energy density rather than a local minimum. To ensure that we retain 
only the physical roots, we have conducted tandem numerical solutions of Eq. H6] and direct 
minimization of the free energy. In the former, instead of solving for a using different 
values of c p , it proves easier to solve Eq. H6] for c p using different values of a. In second 
set of calculations, we have carried out numerical minimization (using Brent's method^) 
of the free energy with respect to a for different values of c p . Results of these two sets of 
calculations are presented in Fig. El Fig. |3k presents the results of the calculations without 
counterions and Fig. |3b corresponds to results with counterions. Solid and dashed lines 
correspond to the roots of Eq. H6] and the numerical minimization, respectively. 

In both the figures, the results of the two sets of calculations match perfectly well except 
in the transition regime, where there are multiple values of the order parameter 5* for a given 
value of c p . In fact, this is the metastable regime for the isotropic-nematic phase transition 
and such a diagram has already been mapped out for the neutral rods using bifurcation 
analysis^. The jump in 5* for the numerical minimization calculations (i.e., dashed lines) is 
a characteristic of the first order isotropic-nematic phase transition. From both the figures, 
it is clear that the isotropic-nematic phase transition takes place at lower monomer densities 
as the linear charge densities of the rodlike polyelectrolytes is increased. Also, in comparing 
the two figures, it is clear that the isotropic-nematic phase transition takes place at higher 
monomer densities in the presence of counterions. These results are consistent with the 
stability analysis of the weakly ordered nematic phase as presented in the previous section. 

We note that the isotropic-nematic transition in this solvated system is actually spanned 
by a region of two-phase coexistence in which a diluted (in polymer) isotropic phase coexists 
with an enriched nematic phase. Such two-phase regions are confined to within the regions 
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of hysteresis shown in Fig. EJ but we have not taken the trouble to elaborate them in this 
figure. The two-phase regions will be revealed when we consider the full phase diagram. 

However, before getting to the full phase diagram, an important remark on the theoretical 
treatment considered in this paper is appropriate. During the computation of the full phase 
diagram, it is found that for some set of parameters a completely ordered nematic phase 
(i.e., the phase for which a — > oo and hence, S = 1) becomes one of the competing phases. 
Physically, the completely ordered phase corresponds to perfectly aligned rods. In the next 
section, we show that Onsager's approach (which is extended to polyelectrolyte mixtures in 
this work) is not able to describe such a coexistence between the completely ordered nematic 
phase and any other phase. 

F. Completely Ordered Nematic Phase for Symmetric Mixtures 



For the completely ordered nematic phase, Onsager's variational parameter, a, diverges. 
In order to evaluate the free energy in this limit, we rewrite Eqs. |4"2T |4"31 using asymptotic 
expansions for the Bessel and hyperbolic functions 34 leading to the result 
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In the limit of a — > oo, the entropic contribution to the free energy of the completely ordered 
nematic phase diverges logarithmically (i.e., /° n = f£ n {a — > oo} — > oo) and the excluded 
volume contribution vanishes (i.e., f° = /^ {a — > oo} — > 0). Furthermore, the limiting 
expression for Eq. 1451 becomes 
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This integral can be evaluated numerically and is found to be negative. Note that Eq. [50] 
is the same as the electrostatic contribution to the free energy considered in Ref.— in the 
context of phase separation of charged aligned needles. 
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The divergence of f° n complicates the evaluation of coexisting densities between the 
isotropic and the completely ordered phases. Also, the vanishing of /£ for any arbitrary 
polymer density highlights the limitation of the Onsager's virial approach to describe phase 
separation when one of the competing phases is fully ordered. 

A more suitable description of such highly ordered phases is presented in Ref.— , where the 
phase separation of a completely ordered parent phase into two completely ordered daughter 
phases having different polymer densities is considered. The phase separation calculations 
were performed using the Percus-Yevick equation of state for hard cylinders. This particular 
description takes into account the higher order terms in describing the excluded volume 
interactions and goes beyond the Onsager approach. However, this study did not consider 
the nematic phase as a candidate in the free energy competition. Such an analysis, while 
possibly relevant in an orientationally constrained situation, could produce unphysical results 
if the system can freely choose the orientation and concentration of both parent and daughter 
phases. 

In the following, in order to avoid the above complications associated with strongly 
ordered nematic phases, we limit our results to regimes where the nematic order parameter 
is not fully saturated at S — 1. 

III. RESULTS 

Using the theoretical approach and free energies presented in the previous section, we 
have investigated the phase behavior of symmetric mixtures of oppositely charged rodlike 
polyelectrolytes. In the following, we consider the possibility of macrophase separation in 
the isotropic phase due to attractive electrostatic correlations between oppositely charged 
polyelectrolytes in addition to phase coexistence between isotropic and the nematic phases. 
Unless mentioned, all lengths are normalized by the monomer length I so that L = Nl, 
with iV the number of monomers. In order to keep the theoretical analysis simple and avoid 
the issue of counterion condensation^, we have focused on the "weakly charged" regime 
corresponding to /#// < 1. Also, we restrict attention to salt-free symmetric mixtures 
here and leave the effect of added salt on the phase behavior for the future. In order to 
identify the role played by the counterions in the salt-free symmetric mixtures, we consider 
two comparison systems - one without any counterions from the polyelectrolytes and the 
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other with counterions. As noted above, the counterion free situation could conceptually 
be realized in symmetric mixtures containing polyacids and polybases, where the charge on 
the polyelectrolyte backbones can be controlled by the pH of the solution. While the theory 
is quite general and a wide variety of phenomena can be investigated, we have chosen to 
further limit the parameter space of our study by fixing the rod length and aspect ratio to 
N = 1000 and d/l = 1/50. 

Prior to presenting results on the full phase diagram, it is illustrative to study macrophase 
separation in the isotropic phase, where we set aside the possibility of nematic order. Even 
for such a simpler situation, there are some key questions that need to be answered. For 
example, it is not clear how the counterions are partitioned between different phases and 
where the phase boundaries (spinodals and binodals) are located relative to the analogous 
phase boundaries in mixtures containing flexible polyelectrolytes. In the next subsection, we 
present our results for phase separation in the isotropic phase and address these key issues. 
In a subsequent subsection, we present the full phase diagram by considering the possibility 
of isotropic-nematic phase behavior in addition to the isotropic-isotropic transition. 

A. Isotropic-isotropic transition 

A salt-free symmetric mixture of oppositely charged polyelectrolytes with counterions is 
a system with five components - two types of polymers, two types of counterions and the 
solvent molecules. However, in the simplified theoretical description presented here, the 
solvent molecules are treated implicitly and by restricting attention to symmetric mixtures 
only, the equations obtained by enforcing the equality of the chemical potentials of the 
polymers in the two phases are degenerate for the two types of polymers. The same is true 
for the counterions. As a consequence, we only have to analyze a pseudo two-component 
system, where the polymers and the counterions need to be partitioned among the coexisting 
phases. Using the same set of arguments, salt-free symmetric mixtures without counterions 
can be treated as a pseudo one-component system. 



21 



1. Effect of counterions 

For the quasi-one and two component systems, we have carried out a direct minimization 
of the total free energy density, obtained by appropriate weighting of two isotropic phases 
and the use of the lever rule 59 . This approach is equivalent to equating the chemical po- 
tential of each component in the two phases and equating the pressure of each phase, but 
replaces solving a nonlinear system by the numerically more robust procedure of minimiza- 
tion. Explicitly, for the quasi-one component system, we minimize the total free energy 
density (f total) °f the phase segregated system with respect to the densities in each phase 
(i.e., two dimensional minimization) after writing 

riso . C P ~ °P f iso r n , °P ~ °p f iso r in ( r-,\ 

J total - r U_ r \J IS/ + r \l_ r lJ i C P/' \° l ) 

where f lso {c l p } and f lso {c p 1 } are the free energy densities of isotropic phases I and II, 
respectively, given by Eq. [13j The parameters c l p , c p , and c l p l (c l p < c p < c l p l ) are the dimen- 
sionless number density of the polyelectrolytes in dilute phase I, overall number density, and 
the number density in the concentrated phase II, respectively. A similar equation can be 
written for the quasi-two component system (i.e. the system with counterions) where we 
note that the total number of counterions is related to the total number of monomers by 
n totai _ 2<jn t L, where n t is the total number of polyelectrolyte chains of one type. For this 
system, three dimensional minimizations of the free energy with respect to the monomer 
densities in each phase along with the counterion density in one of the phases have been 
carried out to compute the coexistence curves. The counterion density in the second phase 
is computed using the lever rule for the counterions. All of the multi-dimensional minimiza- 
tions of the free energy density have been carried out using the simplex method 57 . 

In Fig. HI we present the results of these calculations for different total monomer densities. 
Coexisting monomer and counterion densities in the two phases are presented in Figs. Hk and 
Hb, respectively. In these calculations, the linear charge density a for the polyelectrolytes 
is kept the same so that different total monomer number densities correspond to different 
total counterion number densities. Also, for comparison purposes, the coexisting monomer 
densities in the counterion-free system are also presented in Fig. H^ for a total monomer 
density of 0.002, and are denoted by *. Note that in Fig. 2b the counterion densities on the 
right side of the dashed lines correspond to the densities in the phase with higher monomer 
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density. 

It should be noted that these results are in qualitative agreement with similar calcula- 
tions for flexible polyelectrolytes, where the system phase separates into a very low density 
(supernatant) phase and a dense (coacervate) phase. The asymmetric nature of the co- 
existence curve^ 1 ^ is a result of the long-range electrostatic interactions at play in these 
systems coupled with the reduced translational entropy of polymers relative to solvent and 
ions. Beyond these observations, two important results can be inferred from the figures. In 
Fig. Hk, it is clear that the coexistence regime shrinks continuously with the increase in the 
total counterion density. This implies that the counterions suppress the isotropic-isotropic 
phase transition. This is in agreement with the notion that the counterions need to be par- 
titioned among the phases upon macrophase separation, which is entropically unfavorable. 
Furthermore, from Fig. Hb, it is found that indeed, the counterions get partitioned between 
the two phases and the counterion density is only slightly higher in the coacervate phase. 
Note that this result is in qualitative agreement with the theoretical and experimental re- 
sults reported by VootdM, although the treatment of electrostatics in Voorn's theory is much 
more primitive than ours. 

An important remark regarding the phase diagrams presented in this work is due here. 
Typically, in a two-phase region the boundaries of the phase diagram are the same irre- 
spective of the initial concentration; just the relative amounts of the two phases vary. The 
situation in the presence of counterions is very different. A change in the initial concen- 
tration of polyelectrolytes changes the concentration of counterions in the solution and the 
phase boundaries may shift as described in Fig. HI As different concentrations of the poly- 
electrolytes and their counterions correspond to different states of the system, the diagrams 
presenting the coexisting phases for different initial states should be called "state diagrams" 
in general. However, in this work, we ignore this semantics issue and call these diagrams 
binodals (or coexistence curves). 

2. Comparison with flexible polyelectrolytes: counterion free symmetric mixtures 

In order to compare the phase boundaries of symmetric mixtures containing rodlike poly- 
electrolytes to those of flexible polyelectrolyte mixtures, we have computed the free energy 
of mixtures of flexible polyelectrolytes using an analogous random phase approximation to 
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that employed in the rigid case (see Appendix C). Explicitly, the free energy of a flexible 
mixture is given by 
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Here, if is the excluded volume parameter^ and Q° k is the partition function for a nonin- 
teracting Gaussian chain of length L& = N^h, k = +, — (/& being the Kuhn segment length 
and Nk is the number of segments)^). Furthermore, g(x) = 2(e~ x — 1 + x)/x 2 is the Debye 
function. We note that in the case of flexible polyelectrolyte mixtures the final term in 
F£ (involving the integral) is the well-known Edwards' screening contribution to the free 
energy^ and is negative. For the comparison between the rodlike and the flexible systems, 
this contribution will be ignored and its effect captured by using a renormalized excluded 
volume parameter w r instead of the bare excluded volume parameter w. 

Since our focus here is on the effect of chain flexibility on the electrostatic contribution 
to the free energy, we can identify an appropriate value of w r by forcing agreement between 
the excluded volume contributions to the free energy for the rodlike and flexible symmetric 
mixtures (cf. Eqs. [131 and 152]) . The comparison reveals that w r = irdl 2 /2 is the suitable 
choice, which makes all free energy contributions other than electrostatics identical for the 
rodlike and flexible systems. Using this value for the renormalized excluded volume param- 
eter, the spinodals and the binodals for symmetric mixtures of flexible polyelectrolytes can 
be readily computed using the same approach as used for the rodlike system. To avoid any 
complications arising from the presence of counterions in comparing the phase boundaries for 
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symmetric rodlike and flexible mixtures, we have considered a salt-free and counterion-free 
system. 

For this model flexible polyelectrolyte mixture, the spinodal for the isotropic-isotropic 
transition is given by 
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and c p = 2nL 2 d/Q is the corresponding dimensionless monomer number density. Also, as 
before, we have defined Ib = Ib/1, k = kL, k 2 m = L 2 (k^ + k 2 _), and q = qL. The spinodal for 
the isotropic-isotropic phase transition in a corresponding rodlike symmetric polyelectrolyte 
mixture can be obtained from Eqs. (1561 - |5T|) by replacing the Debye function g (q 2 /6N) by 
the function t lso (q) given in Eq. [T71 

The set of parameters that lead to S\ 80 < corresponds to the regime of instability to 
macrophase separation in the isotropic phase. Note that the electrostatic term D e {c p } is 
positive and hence, the electrostatics drives the macrophase separation in flexible as well 
as rodlike polyelectrolyte mixtures. Also, from Eq. [56] it is clear that the translational 
entropy (which appears as unity in the equation) and the excluded volume interactions 
oppose this driving force (assuming w r > for good solvents). Furthermore, we note 
the prefactor of (4:7rlBZ 2 a 2 N) in front of D e {c p }. This implies that an increase in the 
polymer length, Bjerrum length or the linear charge density leads to a strengthening of 
the electrostatic driving force favoring the macrophase separation. However, note that an 
increase in monomer density leads to screening of the electrostatics, appearing through k sm 
in the expression for D e . This screening effect places an upper concentration bound on the 
unstable regime, ultimately stabilizing a single isotropic phase. 

In Fig. \5[ we have compared the spinodal phase boundaries for flexible and rodlike sym- 
metric polyelectrolyte mixtures in a salt-free, counterion-free situation. It is clear that for 
these weakly charged polyelectrolyte systems (characterized by AnlBO~ 2 p < 1) the isotropic- 
isotropic coexistence regime is broader in the case of flexible symmetric mixtures in com- 
parison with rodlike mixtures. As all the other contributions to the free energies of the two 
systems are the same except the electrostatic contributions, it is clear that the electrostatic 
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driving force is stronger in the case of flexible polyelectrolytes than rodlike polyelectrolytes. 
We note that these numerical results are consistent with the simple scaling arguments pre- 
sented in the introduction. 

B. Full Phase Diagram 

From the results presented and discussions so far, it is clear that electrostatics drives 
macrophase separation in the isotropic phase and also stabilizes the nematic phase. So, 
in principle, isotropic-isotropic, isotropic-nematic and/or nematic-nematic phase transitions 
can take place in symmetric mixtures of rodlike polyelectrolytes. Similar to the isotropic- 
isotropic phase calculations considered above, we consider the extra possibilities of nematic- 
isotropic and nematic-nematic coexistence, where the orientational order parameter in each 
phase is determined during the minimization of the total free energy density of the phase 
separated system. 

For the quasi-one component system without counterions (cf. Eq. ISTj) . we minimize the 
total free energy density written as 

/«, = Jh^/R-*'} +3^4/ Of.""}. (<*> 

Cp C p Cp Cp 

where the variational parameters a 1 and a 11 for the phases I and II, respectively, are ob- 
tained by numerically minimizing (using Brent's method^) the free energy density of the 
phase (/ {c p } and / {cp 1 } for the phases I and II, respectively, given by Eqs. BSl- Hlj) & t the 
given monomer densities (c l and c p l for the phases I and II, respectively). A similar approach 
can be taken in quasi two-dimensional systems that include counterions. The motivation 
behind carrying out such calculations is the fact that all the phase transitions that we have 
considered thus far (i.e., isotropic- isotropic and isotropic-nematic) are subsets of these more 
"general" calculations. Unfortunately, such high dimensional optimizations are computa- 
tionally demanding, so we have looked for opportunities to accelerate the construction of 
full phase diagrams. 

In a complementary set of calculations, we have considered isotropic-nematic and 
isotropic-isotropic phase transitions separately. It was found that the results of the "gen- 
eral" calculations described above exactly match the results of our isotropic-nematic and 
the isotropic-isotropic calculations (data not presented here). Through such calculations, 
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we have also established that there are no regions of nematic-nematic coexistence, so the 
various phase diagrams can be mapped out by tracking individual I-I and I-N boundaries. 

The phase diagrams presented in the remaining sections were obtained by the above 
technique and are presented in Figs. [6] and Ul respectively in the absence and presence of 
counterions. Overall, the results of these calculations are in qualitative agreement with the 
picture presented in Fig. [TJ Some important features of the full phase diagram are discussed 
below. 

1. Phase diagram without counterions 

In Fig. El we have plotted the phase diagram for symmetric mixtures of oppositely charged 
rodlike polyelectrolytes in the absence of counterions. In the figure, the linear charge density 
is varied to explore the effect of electrostatics on the phase boundaries. From Fig. [6^, it is 
clear that the coexisting densities in the two phases decrease from the uncharged case with 
an increase in the charge densities (compare the results for a = 0, 0.02 and 0.04). However, 
a further increase in the charge density leads to three distinct regimes corresponding to 
the low, intermediate and high values of Ib- For very low values of the Bjerrum length 
{Ib/1 —*• 0), the coexisting phases are isotropic and nematic phases. Note that in this regime 
the coexisting densities for the isotropic and nematic phases for the polyelectrolyte systems 
are close to that for the neutral system. 

For high values of the Bjerrum length (close to unity) and high charge densities such as 
a = 0.08,0.1 in Fig. |HK completely ordered nematic phase becomes one of the coexisting 
phases. However, we haven't been able to compute the densities of the coexisting phases in 
this regime due to the numerical issues discussed in subsection F above. 

For intermediate values of Ib, there are regimes (e.g., Ib/1 € 0.1 — 0.6 for a = 0.1 
in Fig. [6k), where isotropic- isotropic coexistence and isotropic-nematic phase separation 
can be separately realized by varying the concentration of the rodlike polyelectrolytes in 
solution. With an increase in Ib in this regime, three phase coexistence (isotropic-isotropic- 
completely ordered nematic) can be realized in these systems. Moreover, the value of Ib 
at which the three phase coexistence takes place is dependent on the linear charge density 
of the polyelectrolytes. In fact, the Bjerrum length at which the three phase coexistence 
takes place decreases with an increase in the linear charge density of the rods (compare the 
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results for a = 0.08 and a = 0.1 in Fig. Ek). Note that in this intermediate regime, the 
entropic effects driving the isotropic-nematic transition are comparable in strength with the 
energetic effects (coming from electrostatics), but the system quickly evolves into a broad 
isotropic-nematic phase envelope upon increasing Is- This leads to a sharp increase in the 
density of the coexisting nematic phase (see the plots for a = 0.08 and o = 0.1 in Fig. [6^). 
A further increase in Ib leads to the completely ordered nematic phase as a coexisting phase. 
In order to keep track of the degree of alignment of the rods in the nematic phase in Fig. [6K 
we have plotted the order parameter in Fig. [6b. 

From Fig. [6b, it is clear that the order parameter increases with an increase in Ib/1, which 
is in agreement with the stability analysis carried out in this paper. Hence, the numerical 
results support the prediction that electrostatics favor orientational ordering. 

2. Effect of counterions 

In contrast to the phase diagram obtained in the absence of the counterions, the phase 
diagram in the presence of counterions depends on the total number density of the rodlike 
polyelectrolytes and correspondingly on the number of counterions in the system. In order 
to conduct a systematic study, we have carried out two sets of calculations. In the first 
set, we vary the linear charge density of the polyelectrolytes while keeping the total number 
density of the rodlike polyelectrolytes fixed at a particular value. In the second set, the total 
number density of the rods is varied keeping the linear charge density fixed at a particular 
value. So, in both the sets, the total number density of counterions is varied. 

Fig. [7J presents the results of the first set of calculations. Comparing Fig. [7J with Fig. [6J 
we observe that the qualitative features of the phase diagram in the presence of counterions 
remain the same as in the absence of counterions. Furthermore, Fig. [8] presents the number 
densities of the counterions in the coexisting phases. It is found that counterions are uni- 
formly distributed between the two phases for these values of linear charge densities and the 
Bjerrum's length. 

The results for the second set of calculations are shown in Figs. M and [TUl where the total 
number density of polyelectrolytes is changed for a particular value of the linear charge den- 
sity. These results show the isotropic-isotropic and isotropic-nematic coexistence regimes 
just like the ones seen in the absence of counterions in Fig. [61 In these calculations, the 
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change in the total number densities of the polyelectrolytes leads to a change in the total 
number density of the counterions. The effect of varying the total number density of the 
polyelectrolytes and counterions on the phase boundaries can be explained as follows. As 
already discussed, the origin of the isotropic-nematic transition for neutral rods near the 
overlap concentration or for very low values of Ib {Ib/1 — > 0) is entropic and is nearly inde- 
pendent of electrostatics. Hence, the presence of counterions does not affect these bound- 
aries (compare the phase boundaries for 2n t L 2 d/Q = 4 and 4.7). However, the origin of 
the isotropic-isotropic transition at low number density of polyelectrolytes is electrostati- 
cally driven and is strongly dependent on the presence of the counterions. The width of 
this particular coexistence regime can be tuned by changing the number of counterions. In 
particular, the coexistence regimes shrink with an increase in the number of counterions 
due to the screening of electrostatic interactions(compare 2n t L 2 d/Q = 0.002,0.01 and 0.1 
in Fig. [H]) • Note that the results for the isotropic-isotropic phase transition in Fig. are the 
same as in Fig. HI 

Figure [10] shows the counterion distribution in the coexisting phases. Consistent with our 
prior results, it is evident that the counterions also phase segregate at high Bjerrum length 
with a higher density in the concentrated (coacervate) phase. 

IV. CONCLUSIONS 

We have studied the phase behavior of salt-free symmetric mixtures of oppositely charged 
rodlike polyelectrolytes using the random phase approximation. In this work, we have 
focused on weak polyelectrolytes in the regime Ib/1 < 1 to avoid complications arising from 
possible counterion condensation. For a variety of symmetric mixtures, we have computed 
the phase boundaries for regions of isotropic-isotropic and isotropic-nematic coexistence. 
We were not able to identify any regions of nematic-nematic coexistence in these symmetric 
systems. 

Our stability analysis and numerical results for coexistence curves reveal that electrostatic 
interactions favor nematic ordering of the rodlike components in solution. Nonetheless, the 
screening of these electrostatic interactions by higher concentrations of counterions weakens 
or destroys this ordering. It is shown that the phase boundaries for symmetric mixtures 
containing oppositely charged rodlike polyelectrolytes are dependent on the electrostatic in- 
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teraction strengths characterized by the linear charge density, a, and the Bjerrum length, 
Ib- In particular, it is demonstrated that at low electrostatic interaction strengths (i.e., 
Ib/1 — > 0, a — > 0), the densities in the coexisting isotropic and nematic phases decrease with 
an increase in the linear charge density of the polyelectrolytes in the absence of counteri- 
ons. However, an increase in the electrostatic interaction strength by increasing a leads to 
three distinct regimes characterized by Ib/1- At low Ib/1 (close to zero), a narrow region of 
isotropic-nematic coexistence prevails, whose origin lies in the entropy of the system. On the 
other hand, at relatively high Ib/1 (close to unity), the completely ordered nematic phase 
becomes one of the coexisting phases. However, its origin lies in the electrostatic attraction 
between oppositely charged polyelectrolytes. At intermediate values of Ib, isotropic-isotropic 
or isotropic-nematic coexistence can prevail depending on the concentration of the polyelec- 
trolytes. Also, in this regime, at a particular value of Ib, three phases (low density isotropic, 
moderate density isotropic, and high density nematic) coexist with each other. The value of 
Ib at which the three phase coexistence takes place depends sensitively on the linear charge 
density of the polyelectrolytes. 

We have also investigated the effect of counterions on the phase coexistence boundaries 
(isotropic-isotropic and isotropic-nematic). Comparison of the results for the systems with 
and without counterions reveals that sufficiently high concentrations of counterions suppress 
both isotropic-isotropic and isotropic-nematic coexistence. A key prediction of theory is the 
result that the concentration of counterions in the dense (or coacervate) phase is slightly 
higher than in the dilute (or supernatant) phase. Also, comparison of the phase separation 
boundaries between comparable rodlike and flexible polyelectrolyte mixtures reveals that the 
isotropic-isotropic macrophase separation regime is broader in the case of weakly charged 
flexible polyelectrolytes. 

Furthermore, in this work, we have limited ourselves to the phase separation regime 
in the mixtures containing weakly charged polyelectrolytes close to the critical point. We 
have found that the critical points for the isotropic-isotropic and isotropic-nematic phase 
transitions exist at very low number density of rods (n/Q ~ 1/L 2 d) and weak electrostatic 
interaction strengths (i.e., Ib/1 <C 1)- At this point, we remark on the range of validity of 
the theory to describe the phase boundaries and some of the future directions. There are 
three key issues, which limits the validity of the theory. In future, we'll extend the theory 
to higher electrostatic interaction strengths by addressing the issues mentioned below. 
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First issue is the use of the random phase approximation (RPA) to compute the con- 
tribution to the free energy coming from the electrostatic correlations. A typical way to 
judge the validity of the RPA is to compare the mean field contribution with the corre- 
lation term in the free energy. Due to the fact that the mean field contribution to the 
free energy coming from the electrostatics is zero and correlation terms beyond the RPA 
are not available, the range of validity of RPA can't be inferred directly. However, general 
consensus 2 ^! 2 ^ - 1 ^ is that the RPA for flexible polymers is valid for concentrations above the 
overlap concentration. In this work, the coacervate phase has concentration above the over- 
lap concentration and the supernatant phase is very dilute. Indeed, the supernatant phase 
is not well described by the RPA. On the other hand, the RPA is suitable for a very good 
quantitative description of the phase boundary describing the coacervate phase. Remark- 
ably, the RPA predicts a supernatant phase with almost zero density in qualitative agreement 
with the experiments^^"— . Furthermore, the RPA describing the coacervate phase bound- 
ary in the case of flexible polyelectrolyte mixtures has been compared extensively with the 
experiments^ and simulations 26 . Indeed, agreement between the theory, simulations and 
experiments is remarkable. In order to go beyond the RPA, Hartree approximation 62 can 
be used, which can provide a more quantiative description for the supernatant phase also. 
However, we leave it for future work. 

Second issue arises due to the use of the Onsager second virial approach to describe the 
steric effects in the case of long rods. It is well-known^ 6 -^ that the approach only works in 
the limit of very long aspect ratio of the rods. In fact, the limit of validity of the Onsager 
corresponds to L/d ^> 10. Furthermore, the second virial approach is strictly valid for low 
number density of the rods. For dense systems, higher order terms needs to be considered. 
For very long rods, the isotropic-nematic and isotropic-isotropic phase transitions occur at 
low enough number densities of the rods (which is of the order of overlap concentration) . This 
is the reason the second virial approach is able to correctly predict these phase transitions at 
low electrostatic interaction strengths. However, an increase in the electrostatic interaction 
strength causes the density of the coacervate phase to increase and the virial approach breaks 
down. This deficiency of the theory can be removed by considering the effect of higher order 
terms 65,66 , which we'll consider in future. 

Third issue is the ignorance of charge renormalization while describing phase separa- 
tion. In this work, the phase separation regime corresponding to very low electrostatic 
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interaction strengths is described assuming that there is no counterion adsorption on the 
polyelectrolytes. In this regime, the issue of charge renormalization during phase separation 
can be safely ignored. However, with the increase in the electrostatic interaction strength, 
counterions can adsorb on the polyelectrolytes and modulate their charge. For this part of 
the phase diagram, the degree of ionization for each polyelectrolyte in each phase needs to 
be computed, while equating the chemical potential and pressure of each component in the 
two phases 6 ^ 1 ^. In other words, the possibility of different degree of ionizations in daughter 
phases needs to be explored. In this work, we have limited the theoretical study to lg/l < 1, 
where issue of counterion adsorption can be safely ignored. 

A related issue is the consideration of ionic clusters (such as dipoles, quadrapoles etc.) 
formed as a result of strong electrostatic attraction between different oppositely charged 
components. It has been shown^ -^ in the literature that the issue of ion-pairing has to be 
taken into account in the case of hard sphere model for simple electrolytes (also known as 
restricted primitive model (RPM)) to correctly match the simulation results in the critical 
regime. This is a manifestation of the fact that the critical point for the RPM, as predicted 
by the Debye-Hiickel theory, corresponds 60 to q 2 l B /a = 16, where q is the charge per sphere, 
l B is the Bjerrum length at the critical point and a is the diameter of the hard spheres. 
For symmetrical electrolytes, a is the same for both kinds of charged spheres. Also, for 
monovalent electrolytes, q = 1 and the critical point exists at very low temperature such 
that l B /a = 16. At such low temperatures, indeed one has to extend the Debye-Hiickel 
theory (which is a RPA like) by including atleast dipolar interactions in the regime near 
the critical point. However, in contrast to the symmetric electrolytic mixtures, the phase 
separation regime close to the critical point in the case of mixtures containing long flexible or 
rodlike polyelectrolytes can be well described within RPA without any consideration of ion- 
ion, ion-polyelectrolyte (which is the same as the counterion adsorption) or polyelectrolyte- 
polyelectrolyte pairs. This is an outcome of the fact that now the phase separation takes 
place at relatively higher temperatures and low densities. In this regime, the mixtures of 
oppositely charged polyelectrolytes behave as weakly correlated liquids 69 . However, if we 
increase the electrostatic interaction strength and go far from the critical regime, there is a 
competition for pairing between different oppositely charged species. A rough estimate for 
the minimum driving force for pairing between the ions is given by Bjerrum's theory^ -^- of 
ion-pairs. According to the theory, life time of the paired state in the case of two oppositely 
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charged ions (while undergoing thermal motion) significantly increases when g 2 /s/r > 2, r 
being the distance between the ions. In this regime, the electrostatic attraction takes over 
the thermal energy of ions and ion-ion pair (or dipole) needs to be considered as a new 
species. This relation can be cast in terms of number density of ions (= n c /Q) by using 
r ~ (n c /f2) -1 / 3 so that q 2 ls/l > l/^e/^) 1 / 3 - Issue of ion pairs/clusters formed as a result 
of binding, which involves charged monomers is more suble compared to the issue of pairing 
in small ions. 

However, we can estimate the regime, where one has to explicitly consider the binding 
between the polyelectrolytes and oppositely charged counterions or polyelectrolytes. Carry- 
ing out a single chain analysis, it has been shown that the binding of counterions on flexible 
and rodlike polyelectrolytes becomes importanti&22 roughly around Ib/1 > 1- Similarly, an 
analysis of the system containing two oppositely charged flexible polyelectrolytes 14 reveals 
that the complexation takes place only for strong electrostatic interaction strengths (i.e., 
<j 2 Ib > 1). For rodlike polyelectrolytes, the electrostatic interaction strength required for 
complexation needs to be stronger compared to the flexible polyelectrolytes due to weaker 
(logarithmic) electrostatic potential for rodlike polymers. These single chain analyses provide 
a clear picture about the dilute solution regime. In the regime above the overlap concen- 
tration for polymers, situation is more complicated and one has to consider the multi-chain 
effects. However, one can carry out a mean-field analysis 2 ^ 3 - to estimate the fraction of 
charged monomers involved in binding (say, T). It can be readily shown^ 3 - that the fraction 
is given by the relation T/(l — T) 2 = (cr/l) [c p /(2Ld)} exp(|.E//csT|), E/ksT ~ Ib being the 
energy gain per pair. In this work, we have considered weakly charged polyelectrolyte solu- 
tions at very low electrostatic interaction strengths so that the fraction of charged monomers 
involved in pairing is close to zero. 

At present we are not aware of experimental data sets sufficient for a comprehensive 
test of the theoretical predictions made here. With this aim, we welcome interactions with 
experimental groups to define appropriate systems and experimental protocols. 
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APPENDIX A : Partition function for mixtures of rodlike polyelectrolytes 



Here, we present the details of our derivation of the free energy for the mixtures of 
rodlike polyelectrolytes in the presence of counterions as described in Section II. Different 
parameters representing the number of rods and counterions, length and diameter of the 
rods, and charge along the rods have already been described in Section II. In terms of these 
parameters, the partition function can be written as 
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Here, Pki^j, n) for k = +, — is the probability of finding a rod of type k oriented along unit 
vector Uj when the director is taken to be along the unit vector n. Subscript Pfe(u, n) under 
the integral symbol means that the integration need to be carried out under the constraint 
that the orientational distribution function for the positive and negative rods are p+(u, n) 
and p_(u, n), respectively. Also, hamiltonian H int is divided into contributions coming 
from the short range excluded volume interactions (H w ) and the long range electrostatic 
interactions (H e ). Furthermore, in the expression for H w , W(rj, Uj, r^, u^) characterizes the 
excluded volume interactions between two rods whose centers are located at rj and r^ with 
their axes oriented along Uj and U&, respectively. Also, l B is Bjerrum length and p e in the 
expression for H e is the local charge density defined as 
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where Z 7 is the valency (with sign) of the charged species of type 7 and Nj, is the number 
of monomers for the rod of type k so that L^ = N^l^, h being the monomer length. Also, 
Sj in the expression for pfc(r) is the contour variable used to locate any monomer along the 
jth rQ( j Q f type k. In order to proceed further, we rewrite Z as 
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where Z is the partition function for a non-interacting system with the same orientational 
distributions (i.e., p + ,p_) and also the normalization factor in the expression for Z, given 
by 
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Using Stirling's approximation for factorials n! ~ n(lnn — 1) and considering different ways 
of distributing n k rods for a given orientational probability distribution function pk(uj,h) 
along the surface of a unit sphere 43 
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Using the Hubbard-Stratonovich transformation^ for the electrostatic part i.e., H e and using 

Fourier transform defined for any arbitrary function /(r) by 
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So far the partition function is exact. However, evaluation of the exact partition function is 
a tedious task. Useful insights can be obtained by invoking the following approximations. 
From Eq. IA-151 
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Note here that the linear terms arising from the expansion of the exponential in Qy vanish 
in the sum due to the presence of oppositely charged counterions, i.e. Z c+ = — Z c _. Fur- 
thermore, approximating the logarithm by the first term in the expansion in Eq. IA-161 is 
strictly valid for dilute concentrations of the counterions so that ny/fl is small. Using this 
approximation, Eq. IA-131 becomes 
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where we have defined k 2 = AnlsYlj' Z^riji/Q so that k _1 is the Debye length. Now the 
integrals over the spatial and orientational degrees of freedom under the constraint of the 
given distribution functions can be carried out using the approximation described in Ref.— . 
To be explicit, 



h {%} 



llfc=+,- llj=l 


f , .s dr,du, 

_ J Pk(Uj,n) J J _ 


exp 


kgT + ^ Z^fc = 4 


,-/ (2^<i Z k (J kPK- q \ 




iw i 


[TS, 


f , ^dv-jdu- 

Jp k (uj,n) J J_ 





37 



llfc=+,- lll = l Jx 



/pfc(Uj-,n) J J 



1 — < 1 — exp 



Hu 



kuT 



exp 



<(■'• 



Efc / T^^qZk^hPk, 



k R T 



. Uk=+,- n"=i [4( Uj -,n) dr J du j] {i - ^p ( 

n fe =+,-nj=i [4( U ,,n) *i du i 
_iiL= 

2 



n"=l Jp fc ( Ui ,n) °M U J / (C^^ 9 Efc ^fc^fcPfc, 



n 2 



+ • • • (A-18) 

n fc =+,-n™=i[4 K ,n) rfr i rfu i_ 

where we have ignored a cross-term between the excluded volume interactions and the elec- 
trostatics part to keep the calculations analytically tractable. Now, assuming that the ex- 
cluded volume interactions between the rods occur independently of each other, second term 
in the series in Eq. IA-18l can be evaluated. 43 Third term in the series is also straigthforward 
to evaluate after plugging the expression for p ksq given in Eq. IA-141 Now, exponentiating 
the series 
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Also, the primed superscript means that the j = k terms are omitted from the double sum. 
By carrying out the Gaussian functional integrals over ip q and subtracting out the free energy 
in the low density limit, Eq. [IT] is obtained. 



APPENDIX B : Limiting cases for the electrostatic part of the free energy of the 

isotropic phase 



Here, we provide a derivation of the approximate form of the electrostatic contribution 
to the free energy of the isotropic phase in the limiting cases of L k — > 00 (cf. Eqs. [22] and 
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. In particular, we focus on the logarithmic corrections to the free energy of the mixture 
of charged rods as already described using scaling analysis in section III Al The electrostatic 
component of the free energy is given by 
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For weakly charged rods so that K^Cl 
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which is Eq. [22]. On the other hand, for strongly charged rods in weakly screened solutions 
k 2 3> 1 and ft — > 0. In this limit, 
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which is Eq. 



APPENDIX C : Partition function for mixtures of flexible polyelectrolytes 



Here, we present the derivation of the free energy for mixtures of oppositely charged 
flexible polyelectrolytes in the presence of their counterions. The flexible polyelectrolytes 
are represented by continuous curves so that Hjk(sj) is the position vector of the s th monomer 
along the j th chain of type k. Furthermore, the contour lengths of chains of type k is taken 
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to be Lk = Nklk, 1^ being the Kuhn's segment length. Accounting for the conformational 
degrees of freedom of the flexible chains by the path integral representation, the partition 
function for this system can be written as 
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Note here that in contrast to the rodlike polymers, the Hamiltonian for the flexible polyelec- 
trolytes has an additional contribution coming from the chain connectivity (Hq) in addition 
to the contributions coming from short range excluded volume (H w ) and long range elec- 
trostatic interactions (H e ). Furthermore, in the spirit of polymer field theories^, we have 
replaced H w by the delta functional form and defined a excluded volume parameter w to take 
care of the short range interactions. Also, the microscopic densities in the above equations 
are defined as 



k B T k B T 

||rfr[p + (r) + p_(r)] 2 
,/5 e (r)p e (r') 



— / dv I dv 

2 



p e (r) = Z + a +P+ (v) + Z_a_p_(v) + ^ %Pj'( r ) 



j'=c+,c- 



■"■k 



3=1 



Nk 



p k (v) = E / ds J S i r ~' R Jk(sj)} } for k = +, 



Pi 



' r 



E*< 



r-ry,j, 



for j' = c+, c — . 



To proceed further, we rewrite the partition function as 
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where Z is the partition function of a mixture of non-interacting chains and counterions. 
Explicitly, it is given by 
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Using the Stirling's approximation n! ~ nlnn — n, we get 
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where _F is the Helmholtz free energy of the mixtures of non-interacting chains and counte- 

n 



rions. Also, Q° k is the partition function of a single Gaussian chain of type k, given by 
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In writing the above equation, we have dropped the index representing the chain number in 
the path integral. Furthermore, 
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where k = +, — and j' = c+, c— . Using Hubbard-Stratonovich transformation^ for the 
excluded volume and electrostatic parts in the Hamiltonian, and using three dimensional 
Fourier transforms as defined in Appendix A, the partition function can be written as 
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where 
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and Qj/ is the partition function for a single small ion of type j', given by Eq. IA-151 Also, 
Qk is the partition function for a single chain of type k, given by 
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Here, pk <q is the Fourier component of the microscopic density of a single chain of type k, 
given by 
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Expanding in powers of pk >q and using translational invariance 
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where g(x) is Debye function, given by 
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Consider the special case of symmetric mixtures so that Z + cr + n + N+g (q 2 N + l 2 + /Q) = 
— Z_a_n_N^_g (q 2 N_l 2 _/Q\. For this particular case, the cross terms containing (pip vanishes. 
Using the approximation for the partition function of small ions as given by Eq. IA-16| the 
partition function for the symmetric mixtures of oppositely charged flexible polyelectrolytes 
becomes 
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where k 2 = 4tiIb X/,-/ =c+ c _ Z 2 ,rij>/Q and K 2 . = AT[lBZla k Nln k /Vt are the same as defined in 
Appendix A for rodlike polyelectrolytes. Carrying out the Gaussian functional integrals and 
subtracting out the free energy in the low density limit, Eq. [52] is obtained. 
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FIG. 1: A qualitative sketch of the phase diagram in a symmetric mixture of oppositely charged 
rodlike polyelectrolytes. The symbols "I" and "N" denote regions of two-phase coexistence involv- 
ing isotropic and nematic phases, respectively. 2n/f2, L and d are the number density of positive 
and negative polyelectrolytes, the length, and the diameter of the charged rods, respectively. The 
y-axis is the ratio of the Bjerrum length lg to the monomer length I and is a measure of electrostatic 
strength that is inversely proportional to temperature. 
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(b) 
FIG. 2: Phase boundaries (spinodals) for the stability of a weakly ordered nematic phase for a 
symmetric mixture of rodlike polyelectrolytes in the absence (a) and the presence of the counterions 
(b). On the right hand side of each boundary, the nematic phase (denoted by "N") is stable in 
comparison with the isotropic phase. These results are obtained for N = 1000, d/l = 1/50. 
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(b) 
FIG. 3: Effect of the linear charge density a on the isotropic- nematic phase transition (charac- 
terized by jump in order parameter) in salt-free symmetric mixtures of oppositely charged rodlike 
polyelectrolytes. Figs, (a) and (b) correspond to trends in the nematic order parameter S with 
respect to dimensionless rod concentration 2nL 2 d/0, for symmetric mixtures without and with 
counterions, respectively. In each figure, the solid lines present the solutions obtained by solving 
the nonlinear equations for the optimal variational parameter in Onsager's trial function and the 
dashed lines are the results of numerical minimization of the free energy with respect to the same 
variational parameter. The curves correspond to N = 1000, d/l = 1/50 and lg/l = 0.7. 
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(b) 
FIG. 4: Coexistence curves (binodals) for isotropic-isotropic phase separation in salt-free symmetric 
mixtures of oppositely charged rodlike polyelectrolytes. Polyelectrolyte density and counterion 
density in the coexisting phases are presented in Figs, (a) and (b), respectively. 2n t /VL corresponds 
to the total number density of rodlike polyelectrolytes and these curves were obtained using N = 
1000, d = 1/50, a = 0.1. In Fig. (a), the plot with 2n t L 2 d/n = 0.002* represents the binodal 
for a symmetric mixture without any counterions for comparison purposes. Also, in Fig. (b), the 
counterion density in the dense phase corresponds to the right hand side of the dashed lines. In Figs, 
(a) and (b), dashed lines indicate the total number densities of polyelectrolytes and counterions, 
respectively. 




2nL 2 d/Q 
FIG. 5: Comparison of the envelopes for isotropic-isotropic phase separation in symmetric mix- 
tures of the oppositely charged flexible and rodlike polyelectrolytes. Dashed curves correspond to 
the spinodals and solid curves represent the binodals. The curves correspond to N = 1000, d = 
1/50, a = 0.1. "I" and "I-I" denote one-phase isotropic and two-phase isotropic-isotropic coexis- 
tence regimes, respectively. 
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(b) 
FIG. 6: Phase diagram for symmetric mixtures of oppositely charged rodlike polyelectrolytes in 
the absence of counterions and for varying linear charge density a. Figs, (a) and (b) correspond 
to the coexisting number densities of rods in the two phases and the nematic order parameter 
S, respectively. The symbols "I" and "N" denote the isotropic and nematic phases, respectively. 
These curves correspond to N = 1000, d/l = 1/50. Dashed lines in Fig. (a) correspond to the 
boundary beyond which the isotropic phase becomes unstable and the completely ordered nematic 
phase becomes one of the coexisting phases. 
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(b) 
FIG. 7: Phase diagram for symmetric mixtures of oppositely charged rodlike polyelectrolytes in 
the presence of counterions, at fixed total polymer density, and for varying linear charge density 
a. Figs, (a) and (b) correspond to the coexisting number densities of the rods in the two phases 
and the nematic order parameter S, respectively. Symbols "I" and "N" denote the isotropic and 
nematic phases, respectively. These curves correspond to N = 1000, d/l = 1/50, 2ntL 2 d/Q = 4.7. 
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FIG. 8: Counterion densities in the two coexisting phases shown in Fig. [7J The counterions are 
uniformly distributed between the two phases for these set of parameters. 




I + I 



0.25 
0.2 

0.15 
0.1 







I + N 



6 8 

2nL 2 d/Q 



■2n t L 2 d/Q = 0.002 
2n t L 2 d/tt = 0.01 
■2n t L 2 d/tt = 0.1 
■2n t L 2 d/Q = 4.0 
■2n t L 2 d/Q = 4.7 




FIG. 9: Effect of the total number density of rodlike polyelectrolytes on the phase diagram in 
symmetric mixtures in the presence of counterions and at fixed linear charge density. These plots 
correspond to N = 1000, d/l = 1/50 and a = 0.1. Dashed lines in Fig. (a) correspond to the 
boundary beyond which the isotropic phase becomes unstable and the completely ordered nematic 
phase becomes one of the coexisting phases. 
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FIG. 10: Counterion densities in the coexisting phases corresponding to I-I and I-N phase transi- 
tions as shown in Fig. In these figures, the right branch of each coexisting curve corresponds 
to the counterion density in the phase dense in rodlike polyelectrolytes (coacervate) . It is shown 
that for strong enough electrostatics, the counterions also phase segregate and the counterion den- 
sity is higher in the dense phase (see the results for 2ritL 2 d/£l = 0.1). These plots correspond to 
N = 1000, d/l = 1/50 and a = 0.1. 
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